Transition Path, Quasi-potential Energy Landscape and Stability of Genetic Switches 
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One of the fundamental cellular processes governed by genetic regulatory networks in cells is the 
transition among different states under the intrinsic and extrinsic noise. Based on a two-state genetic 
switching model with positive feedback, we develop a framework to understand the metastability 
in gene expressions. This framework is comprised of identifying the transition path, reconstructing 
the global quasi-potential energy landscape, analyzing the uphill and downhill transition paths, etc. 
It is successfully utilized to investigate the stability of genetic switching models and fluctuation 
properties in different regimes of gene expression with positive feedback. The quasi-potential energy 
landscape, which is the rationalized version of Waddington potential, provides a quantitative tool 
to understand the metastability in more general biological processes with intrinsic noise. 

PACS numbers: 87.18Cf, 02.50.Ey, 82.39.-k, 87.17.Aa 



In a cell, the reactions underlying gene expression in- 
volve small numbers of molecules, such as DNA, mR- 
NAs, and transcription factors, so the stochasticity in 
gene regulation process is inevitable even under constant 
environmental conditions [TJ [2]. Recent progresses in 
single-cell observations and analysis demand new quanti- 
tative methods to characterize the regulatory mechanism 
in both prokaryotes [3] and eukaryotes [4] from the cel- 
lular stochasticity. 

Previous kinetic studies of cellular stochasticity have 
been formulated by using the generating function [5 j , sys- 
tem size expansion [6] [7], large deviation theory [8j [9], 
or by employing WKB approximation to the chemical 
master equations (CMEs) [T0UT2] . etc. However, only 
few of them take account of transcriptional noise explic- 
itly. Some recent studies have shown that correlations 
between mRNA and protein levels do not always per- 
form equally well in revealing genetic regulatory rela- 
tionships [T3l H3], and the consideration of mRNA has 
intensive effect on the switching times p~5j [16]. Since 
Waddington's "epigenetic landscape" proposed in 1957 
[T7] , the energy landscape have been widely used to pro- 
vide pictorial illustration of the dynamics and evolution 
of genetic regulatory systems [2 , 8 , 18 . Thus it is impor- 
tant and natural to develop a general methodology which 
can effectively determine the key features of a gene ex- 
pression system, such as constructing the corresponding 
"Waddington potential" , identifying the transition paths 
between metastable states and computing the transition 
rates, etc. 

In this letter, we present a framework to understand 
the metastability of the genetic switches in gene expres- 
sion based on the large deviation theory for Markov pro- 
cesses [HHZT. By explicitly taking into account the 
mRNA noise, we obtain the most probable transition 
paths for off-to-on and on-to- o/f genetic switches through 
the geometric minimum action method (gMAM) [22] . 



Furthermore, we reconstruct the global quasi-potential 
energy landscape, which is the rationalized version of the 
Waddington potential, via the computation of the local 
quasi-potentials. We analyze the properties of transition 
paths, discuss their relation to the quasi-potential energy 
landscape and the deterministic dynamics, and compare 
the differences between our approach and others in pre- 
vious literatures. Our analytical results agree well with 
the Monte Carlo simulations. We also discuss the choice 
of the system size and its finite effect via the transition 
path theory (TPT) [23]. From the authors' opinion, this 
framework is generally applicable for studying transitions 
between stable-saddle- stable fixed points with jump type 
noise generated by Gillespie's birth-death dynamics [24] . 
It is successfully utilized to investigate the stability of ge- 
netic switching models and fluctuation properties in dif- 
ferent regimes of gene expression with positive feedback, 
which leads to interesting biological insights. 

Let us focus on a two-state gene expression model 
in Fig. [IJa) following the previous study [12], where 
the transitions of active and inactive promoter states 
are controlled by the protein via a positive feedback. 
The transition rates between active and inactive states 
are k on (n) = f(n) and k ff(n) = g(n), where n is 
the protein copy number, on and off denote the high 
and low states of protein and mRNA respectively. Here 
the transition rates are in the form of Hill-function as 
f{n) = k^ in + (k^* - k^ in )n hl /(n£<J + n hl ) and g(n) = 
fc max _ ^max _ k^ n )n h * / (n 1 ^ + n h2 ), where n 50 is the 
curve's midpoint, and hi = /i2 = h are chosen for sim- 
plicity. 

A deterministic description of the genetic switching 
model in Fig. [I] is given by the equations 

M = af(N)/[f(N) + g(N)} - 7 M, N = jbM - N, (1) 

through quasi-steady state approximation, where M,N 
denote the mean number of mRNAs and proteins, re- 
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FIG. 1. (color online), (a) Model for gene-expression switch- 
ing network following Assaf et al [12]. Promoter transi- 
tions are regulated by the feedback functions k on and k Q ff. 
Transcription, translation and decay of mRNA and protein 
are modeled as first-order reactions with rates a, 76, 7 and 
1 (rates are rescaled by the protein decay rate), (b) Visit- 
ing frequency distribution of mRNA and protein with a long 
time Monte Carlo simulation (switches occur 2000 times). 
Darkness of points shows the number of visits with suit- 
able smoothing. The blue solid line is the separatrix be- 
tween two basins of attraction of on and off states. In (b), 
K = ab = 2400,6 = 22.5, h = 2,n 50 = 1000, k^ in = kf n = 
a/100, A^ ax = ^ ax = a and 7 = 2. 

spectively. The corresponding stochastic description is 
governed by the CMEs [7] 

Qm,n = -g(jl)Qm,n + /W^m.n + [A + a(E~ 1 - l)]Qm,n- 

where P m ,n denotes the probability distribution function 
(PDF) of inactive DNA/promoter state with m mRNAs 
and n proteins at time t, while Q m , n denotes the active 
DNA/promoter state. Here we employ the notation for 
raising operator W n acting on f(n) as W n f(n) = f(n+j), 
and A = (E* - l)n + 7(E^ - l)m + 7^m(E- 1 - 1) 
is a birth-death operator related to the inactive pro- 
moter. Equations exhibit the bist ability but ignores 
the noise. Under the deterministic description, once the 
system settles in one of its two attractive fixed points, it 
will stay there forever. However, in the presence of in- 
trinsic noise, the system will fluctuate around its attrac- 
tive fixed points and switch between its two metastable 
states on a large timescale (see Fig. [IJb)). In the regime 
of interest (see also [12 ), we fix parameters 7,6 and 
let a = Kb -1 . Here K plays the role of system size 
[7J [19] . We will let K goes to infinity and fix the ratios 
n^/K, k^ in /K, k^/K, kf n /K, kf**/K in the limiting 
process. This choice gives the rationale that in which 
sense the ODEs is the mean field limit of the CMEs 
([2} (cf. the Supplementary Material (SM)). 

The large deviation theory gives a good quantitative 
description of rare events [I9j [22j [25] . The most probable 
transition path is given by minimizing an action func- 



tional characterized by a Lagrangian L. For the Gille- 
spie's birth-death dynamics, L has no closed form and 
only its dual Hamiltonian can be given as 

N 

H(x,p) = J2^)(^ j -^ (3) 

where aj is the rate (or propensity) and i/j is the state- 
change (or stoichiometric) vector, j = 1, 2, . . . , iV. How- 
ever, this description has difficulty when taking DNA into 
consideration, for normally there are too few DNA copies 
in a living cell that the straightforward application of the 
above formulation is not valid. 

Now we adopt Assaf et al's approach [12] by consider- 
ing the stationary distribution where P m , n = Qm,n — 
and eliminating Qm,n to obtain 

= {A + g(n)- Y \K + a(E" 1 - l)][/(n) - A]}P m , n . (4) 

Define the concentration variable x = m/K, y = n/K 
and the quasi-potential S(x, y) [20 , we plug WKB ansatz 

Pm,n = P(x, y) ~ exp[-KS(x, y)] (5) 

into Eq. Q to get a stationary Hamilton- Jacobi equa- 
tion H(x,y,d x S,d y S) = to the leading order with the 
Hamiltonian 

H = A + g(y)- 1 [A + b- 1 (e^ - 1)] [f(y) - A] , (6) 

where the function A — A(x,y,p x ,p y ) = y(e~ Py — 1) + 
jx(e~ Px — 1) + jbx(e Py — 1), and the feedback functions 
are now rescaled as f(y) = f{y)/K,g(y) = g(y)/K. The 
classical Hamilton- Jacobi theory enables one to solve the 
quasi-potential S(x,y) with different formulations such 
as the variational methods, integrating Hamiltonian dy- 
namics or directly solving PDEs, in which we have the 
connection p x = d x S and p y = d y S between the mo- 
menta and the quasi-potential [26]. It is worth asking 
whether the choice of the large parameter K affects the 
final results since any choice is artificial in practice. An 
affirmative answer is given in SM that only the scaling 
matters and the final systems are equivalent with respect 
to different choices of the large parameter K. 

Now we use the powerful gMAM algorithm [22] to com- 
pute the quasi-potential by minimizing the action func- 
tional with the obtained Hamiltonian (|6| (some details 
about the gMAM can be referred to SM) . We choose the 
two stable points solved by Eqs. as our starting and 
ending points, then compute the switching paths start- 
ing from either of the two states (see Fig. [2j. These 
switching paths predict well with the switching trajecto- 
ries obtained from the MC simulations. 

Figure [2] shows clearly that when switch occurs, the 
switching trajectory prefers to be around the most prob- 
able path characterized by the Hamiltonian (|6|. Fur- 
thermore, the o ff-to- on and on-to-off paths are not iden- 
tical, suggesting that the switching paths are irreversible, 
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FIG. 2. (color online). Switching paths (a) from off to on 
state (red solid curve) and (b) from on to off state (purple 
solid curve) and MC simulations for both switching trajecto- 
ries. We take the two stable fixed points in the deterministic 
dynamics as the starting and ending points. Darkness of the 
shading points represents the number of visits for reactive tra- 
jectories with smoothing. The results are obtained from 1000 
independent long time MC simulations. Here, the parameters 
are the same as in Fig. [I] 



which is different from the original Waddington picture 
[T7] . The irreversibility is fundamental in chemical reac- 
tion kinetics due to the non-gradient nature of the dy- 
namical system. Figure [2] also tells that both switch- 
ing paths pass through the same bottleneck, i.e. the 
saddle point given by Eqs. 0. To further character- 
ize the switching path, we note that the transition path 
is also given by the Hamilton's equations x = V p i7, 
p = —V X H, where x = (x,y),p = (p x ,Py)- Based on 
the fact H(x : y, 0,0) = 0, we obtain V X H = when 
p = 0. At the saddle point in any transition path, we 
have p = [22] , and thus p = along the whole downhill 
path. With this result we obtain the downhill equations 
x = V p i^(a?,0), which exactly corresponds to the de- 
terministic dynamics with a renormalization factor 
9(y)/[f(y) + d(y)] °f time si nce this only involves the 
DNA inactive state. This fact explains that after climb- 
ing the saddle point the biological system relaxes to its 
attracting state fast without cost any action. The uphill 
path is then given by x = \7 p H(x, \7 X S) which is in gen- 
eral different from the downhill path unless the system is 
of gradient type. 

Besides giving the minimum action path, gM AM also 
contributes the action value at each point in the path, 
which enables us to calculate the mean switching time 
(MST) r from either stable state. Denote the quasi- 
potential energy barrier AS'on = So — S on , where So is 
the action at the saddle point and S on is the action at 
the on state. Then the MST r on from on-to-off transition 
can be roughly estimated from an asymptotic analysis 
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where T on is a prefactor which is independent of K in 
many cases. Although for one dimensional systems the 
prefactor of MST can be obtained [27], there are no avail- 



able results in high dimensions because of the geome- 
try problem in more than one spatial dimension and the 
non-gradient nature of the system [28] [29] . Therefore we 
compare the MC simulations with the exponential time 
part and adjust the prefactor T on to fit the numerical re- 
sults. Figure [3ja) and (b) demonstrate the MST versus 
the variation of parameters b and 77,50, respectively. It 
shows that the MST is excellently predicted by Eq. ^ 
up to a slowly varying prefactor. The analysis for r ff 
is similar. Our results in Fig. |3ja) show that the MST 
from both off-to-on and on-to-off states decrease rapidly 
when translation rate b is increased. And when n^o in- 
creases in the considered interval, MST from off to on 
state increases exponentially, while MST from on to off 
state decreases exponentially. When n§o is about 1004, 
these two MSTs are equal. These results may provide 
hints to the range of kinetic reaction rates. 
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FIG. 3. (color online). The Mean Switching Time (MST) 
from both states as a function of (a) translation rate b (7 
is hold constant) and (b) Hill-type function curve's midpoint 
7150. The gMAM results with numerical prefactor from off 
to on states (red solid line) and from on to off states (blue 
dashed line), while MC simulations (•) and (o), respectively. 
Prefactor from off to on state were (a) 31.08, (b) 27.61 and 
from on to off states (a) 4.82, (b) 7.94. Here h = 2, K = ab = 
2400, k^ in = kf n = a/100, k^ ax = k? ax = a, 7 = 2. And for 
(a) n 50 = 1000, (b) b = 22.5. 



Based on the obtained quasi-potential values for each 
point from on and off states as starting point, we can re- 
construct the global quasi-potential energy landscape for 
genetic switching model. The two and three dimensional 
view of the quasi-potential energy landscape is shown 
in Fig. [4] and the reconstruction details can be found 
in SM. The reconstruction for more general stochastic 
dynamical systems can be referred to [20] [3D]. In Fig. 
[4j we observe that the on and off states correspond to 
two local minimum on the quasi-potential energy land- 
scape, the saddle of the deterministic dynamical system 
exactly corresponds to the saddle point on the quasi- 
potential energy landscape, too. As discussed in [20], we 
can obtain that the most probable uphill path satisfies 
(f = a ((p) • \7U(cp) + 1 ((f), and the most probable down- 
hill path satisfies ij> = b(r/>) = -a(ip) • VU(ip) + 
if the underlying stochastic dynamics has the form x = 
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FIG. 4. (color online). Quasipotential energy landscape of 
the whole genetic switching system with (a) two and (b) three 
dimensional view as well as switching paths between two sta- 
ble fixed points. Each path passes through the saddle point. 
Here, K = ab = 2400,6 = 22.5, h = 2,n 50 = 1000, k^ in = 
kf n = a/100, £^ ax = k? ax = a and 7 - 2. 



b(x) + y/e<r(x) • w and the drift b has the decomposition 
b(x) = -a(x)-VU(x)+l(x) such that l(x) • VU(x) = 0, 
where w is the standard temporal Gaussian white noise 
and the diffusion matrix a(x) = a(x)-a(x) T . The math- 
ematical derivations and observations in Fig. [4] show 
that the similar picture still holds for the chemical jump 
processes but the argument by simple orthogonal type 
decomposition of the drift is no longer valid if we re- 
call that the uphill dynamics is x = \7 p H(x,\7 x S). In 
general, this form does not permit to specify some ma- 
trix a and make a meaningful decomposition because 
of nonlinearity. But it is instructive to remark that 
x = \7 p H(x,\7 x S) « V p H(x,0) + V 2 pp H(x,0) ■ \7 X S 
when x is close to the critical points with property 
\7 X S = 0. Based on the fact that x = V p H(x,0) gives 
the mean field dynamics, the above derivation actually 
states that the behavior of the system is well approxi- 
mated by the chemical Langevin process when the tran- 
sition path is close to the metastable states and saddle 
points. More detailed discussions on this point may be 
referred to SM. 

The quasi-potential energy landscape not only pro- 
vides the pictorial illustration for the dynamical tran- 
sitions in genetic switching models, it also includes more 
quantitative information to understand the metastabil- 
ity in gene expressions in single-cell observations PQ, 
such as the fluctuation property like the ratio of the 
standard deviation to the mean (SMR). The SMR for 
our genetic switching model can be obtained by the fol- 
lowing calculations. From Eqs. Q and ([5| we know 
that the contribution of the inactive and active promoter 
states to the QSD has similar forms P(x,y),Q(x,y) ~ 
exp[—KS(x,y)]. Thus the whole QSD V{x,y) reads 
V{x,y) = P(x,y) + Q{x,y) ~ e -KS(x,y)^ ^ e can ex _ 
pand S(x,y) in the vicinity of (x on ,y on ) up to second 



order thus get the Gaussian approximation 

V{x, y) ~ (27r) | E| i ex P { " \(x-^~Hx-v) T }- (8) 

Here, x = (x, y), fi = (x on , y on ), E = (6^)2x2, and \^\ is 
the determinant of matrix X. Eq. ^ holds only in the 
vicinity of the on state with standard deviations a m = 
(KS'> y /\?:\)KcTn = (tfS^/lED*. With the a m and a n 
above, we can easily obtain the SMR as shown in Fig. [5j 




Z , ^UUU JWU tuw 

mean number of protein mean numD er of protein 



FIG. 5. (color online). The SMR vs mean number of protein 
induced by varying promoter activation rate k on (a) and tran- 
scription rate a (b). Analytical results of eukaryotic model 
(red solid line) and prokaryotic model (blue dashed line), 
while MC simulations (•) and (o), respectively. 

We distinguish among several mechanisms of 
DN A/promoter activation in genetic switches by 
selecting different promoter transition rates a and 
mRNA degradation rates 7, which is similar to the 
previous work of Raser and O'Shea [3]. The above 
model in Fig. [I] is set with the parameters a=106.7 and 
7 = 2.0. We set a slow chromatin-remodeling eukaryotic 
model with k^ ax = k™ ax = a/10, k^ in = kf m = a/1000, 
and 7 = 1.5, while a prokaryotic model is set as 
fc max = kf ax = 10a, k^ in = k^ in = a/10, with a higher 
mRNA degradation rate 7 = 15. The results in Fig. [5] 
show that the noise level in eukaryotic model is almost 
always bigger than that in prokaryotic model. And the 
relationship between SMR and mean expression level is 
different from the no- feedback models [4 j, especially for 
the changing of transcription rate a. It is interesting to 
find that when increasing the transcription rate a, noise 
level of on state decreases in the model with positive 
feedback, while protein noise increases in one-state 
chromatin-remodeling eukaryotic model without positive 
feedback (Case I in Raser and O'Shea's work). In SM, 
we also compare the SMR curves of different promoter 
transition rates a with fixed 7 = 2.0 (Fig. S4) and 
provide the stochastic features of mRNA (Fig. S3). 

Our methodology above is based on large deviation 
theory, which requires that the system size tends to in- 
finity. For more general and realistic cases, the system 
size is finite, thus the transition path theory (TPT) can 
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FIG. 6. (color online). Switching paths compared with 
gMAM result. gMAM result is colored green in both sub- 
figure, (a) from off to on state and (b) from on to off state 
and MC simulations for both switching trajectories. We take 
the nearest 4 lattice points around starting and ending point 
in gMAM as the starting and ending set. Red line means the 
system is in DNA active state while blue line means DNA is 
inactive. The parameter setup is the same as in Fig. [I] 



be applied to find the most probable transition paths [23 . 
Here we choose the nearest 4 lattice points around start- 
ing and ending point in gMAM as our starting and end- 
ing sets and perform the TPT computations (see Fig. p 
The implementation details can be found in SM. Fig. 
shows that TPT result fits MC simulation slightly better 
than gMAM. TPT result in Fig. |6ja) also provides inter- 
esting details to understand the uphill trajectory from 
off to on state. In the beginning the off-to-on switch, 
the system is mostly at the active DN A/promoter state 
thus the mRNA level increases before the translation of 
protein. While in the beginning of on-to-off switch in 
Fig. [6](b), we observe opposite mechanism with inactive 
DN A/promoter and first decrease of mRNA. This ex- 
plains why the uphill transition path from off to on state 
is convex and it is concave from on- to- off state. 

In this letter, we have presented an analytical frame- 
work to reconstruct the quasi-potential energy landscape 
of genetic switching system while explicitly taking mRNA 
noise into account. We further analyze the properties 
of the transition paths and clarify the relation with the 
gradient-curl decomposition in the previous literatures. 
This global potential, which is a rationalized version of 
Waddington potential, provides a quantitative tool to 
understand the metastability in more general biological 
processes with intrinsic noise. The applications to other 
biological systems such as complex cellular decision mak- 
ing process and the development process of cells will be 
investigated in the future. 
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